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Abstract. We introduce a discrete time microscopic single particle model for kinetic 
transport. The kinetics is modeled by a two-state Markov chain, the transport by 
deterministic advection plus a random space step. The position of the particle after 
n time steps is given by a random sum of space steps, where the size of the sum 
is given by a Markov binomial distribution (MBD). We prove that by letting the 
length of the time steps and the intensity of the switching between states tend to zero 
linearly, we obtain a random variable S(t), which is closely connected to a well known 
(deterministic) PDE reactive transport model from the civil engineering literature. 
Our model explains (via bimodality of the MBD) the double peaking behavior of the 
concentration of the free part of solutes in the PDE model. Moreover, we show for 
instantaneous injection of the solute that the partial densities of the free and adsorbed 
part of the solute at time t do exist, and satisfy the partial differential equations. 
Key words. Markov binomial distribution, reactive transport, kinetic adsorption, solute 
transport, multi-modality, double-peak. 
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1. Introduction 

We consider a mathematical model for the displacement of a solute through a medium 
which apart from a constant flow (advection) and a dispersion (diffusion) interacts with 
the medium by intermittent adsorption (the kinetics). Our goal is to connect a stochastic 
single particle model to the well known deterministic model which describes this process 
by a pair of partial differential equations. 

In Section[2]we give an introduction to the deterministic reactive transport model (as e.g. 
in |10j ) characterized by a pair of partial differential equations. 

In Section|3]we give our simple discrete time microscopic single particle stochastic reactive 
transport model. In Section |4] we calculate the probability generating functions of the 
Markov binomial distribution (MBD) which is described in Section [3J These are helpful 
to consider the convergence of our simple discrete time stochastic model by letting the 
time step go to zero. This will be discussed in Section [5] In Section [6] we compare our 
discrete time model with the obvious continuous time model. 
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In Section [7] we show for instantaneous injection of the solute that the partial probability 
densities of the free and adsorbed parts of the solute do satisfy the PDE's defined in 
Section [2] In Section [5] we compute the means and variances of our stochastic reactive 
transport model. Actually our formula fills a gap in jlOJ : since the authors erroneously 
state that the variances are linear in the initial distribution, they only give the result for 
two initial distributions (this might be connected to their formula (22), which is incorrect). 
In Section Owe study the probability density function of our stochastic reactive transport 
model. This gives us a new and more precise point of view at the double peaking behavior 
in the concentration of the free part of the solute discussed by Michalak and Kitanidis in 

2. The PDE reactive transport model 

We describe shortly the model used by Michalak and Kitanidis in [10] (see [9] for a more 
extensive treatment). Given is a solute that has a sorbed part that does not move, and 
a free part that moves in the x-direction by advection and dispersion. Let C F (t, x) and 
Ca (t, x) denote the concentration functions of the free and the adsorbed part of the solute 
at time t at position x. By applying mass conservation and Fick's law one can set up the 
following pair of differential equations: 



dC F (t,x) dC A {t,x) d 2 C F {t,x) dC F (t,x) 



(1) 



dt 



dt 

dC A (t,x) 
dt 



D 



dx 2 
/iC A (t,x) 



dx 

XC F (t,x). 



Here D is called the dispersion coefficient and v the advection velocity. The parameters 
A and \x denote the rates of changes as described in Figure (TJ with A for the change from 
free to adsorbed and fj. for the change from adsorbed to free. 



D 

- Dispersion - 




Figure 1. The schematic description of the kinetic transport model. 



The initial and boundary conditions are given by 

dC T (t,x) 



C T (0,x) = v r 5(x), lim C T (t,x) — lim 



dx 



= for t > 0, r € {F, A} 
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where (vf, va) is a probability vector and 5 the Dirac delta function. 
Michalak and Kitanidis have a slightly different set up, where the basic quantities are the 
aqueous concentration C and the contaminant mass sorbed per mass of aquifer solids S. 
The connection is given by 

C F = i]C, C A = pS, 
where rj is the porosity and p mass of aquifer solids per total volume. 
Also, Michalak and Kitanidis do not directly use A and (i, but rather consider a distribution 
coefficient and a mass transfer coefficient k, which are given by 

A = — — -k, p = k. 

7] 

The main goal of the authors of [10] is to obtain closed form expressions for the m th 
normalized moments for the free and adsorbed phase, defined by 

+oo 

M^\t) = ^— J x m C T (t,x)dx, re{F,A}, 

— OO 

" — + few* 

These moments (for m — 1 and m = 2) arc obtained in [lOj by taking Fourier transforms 
in the partial differential equations (p}, and differentiating. We copy here the formulaQ 
from ([TD], page 2136) for the normalized second central moment p^it) where the solute 
is in the free phase both at time and at time t: 

t 2 Av 2 /3 (8 - if ( 2D 2v 2 (3 \ 



Unit) = 5 ^ + t 

(2) + tA 



(B + i f (1 + 13 A) 2 + l k (/3 + i f J 

'4v 2 f3 (-B 2 A~ B 2 - (3 + 1) 2D f3(/3-l) 



\ k(l + (3Af((3 + lf (P + 1)(1 + BA)J 
2v 2 (3 (1 - A) (3 B 2 A - 3 - (3 (A © 1)) AD (3 (1 - A) 



k 2 (l + (3Af(B + lf k(l + (3A)(B + lf 

Here Michalak and Kitanidis have made the following abbreviations: 

/?= — = -, A= A(t) = exp(-(/3 + l)kt) = exp(-(A + p)t). 
VP 

3. A SIMPLE STOCHASTIC REACTIVE TRANSPORT MODEL 

We describe the behavior of a single particle in the solute. Time t is discretized by choosing 
some n, and dividing [0,t] into n intervals of the same length 

At = t/n. 

We suppose in such an interval of length At that the particle can only be in one of the two 
states: 'free ' or 'adsorbed ', which we code by the letters F and A. The particle can only 
move when it is 'free ', and in this case its displacement has two components: dispersion 



Hhe © is a + sign in ([10]), but should be a —sign 
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and advection. Let Xk , k > 1 be the displacement of the particle due to the dispersion the 
fcth time that it is 'free'. We model the Xk as independent identically distributed random 
variables satisfying 

(3) E u [X k ] = 0, E„ [X%] = 2D At and E„ [X%] = o(At) as At 1 0, 

where D > 0, and v = (i^,i/a) is the initial distribution describing the state of the 
particle at time 0. When the particle is free during the interval [(k — l)Ai, kAt] for some 
k, the displacement due to advection is given by vAt with v the (deterministic) advection 
velocity. 

In order to model the kinetics, let {Yk,k > 1} be a process taking values in {F, A} (we 
will make a choice for {Yfc} below), and let 

n 

K n = ^2 1 { y K= F } 

k=l 

be the occupation time of the process {Yfe} in state F up to time n. 



X 




At 2At 3A/ lAt 5At ' • ' t-At t 



FIGURE 2. The position S n (t) of the particle at time t = nAt with 

Yi =A,y2 = F ) y3 = F ) y4 = A,y 5 = F,--- ,f„ = f. 

Now let S n (t) be the position of the particle at time t = nAt. Then by the above (see 
also Figure [5]) we can write S n (t) as 

S n (t) = J2{X k + vAt). 

k=l 

Here we assume that K n is independent of the dispersion Xk, k = 1, . . . , K n . 
We want to compare our stochastic model with the PDE-model of Michalak and Kitani- 
dis from Section [5] Since these authors consider the solute with given states ('free' or 
'adsorbed') at time t, we need to consider the conditional random variables S^(t) and 
S^(t), i.e., the position of the particle at time t given that it is 'free' and 'adsorbed' 
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respectively at time t = nAt. Let be the random variable K n conditioned on Y n = r 
with r 6 {F, A}, i.e., counts the number of intervals [(k — l)At, fcAi], 1 < k < n where 
the particle is free, conditioned on the particle being in state r in [t — At,t]. Then S£(t) 
can be written as 

A'; 

S£(t) = £(X fc + t;At). 
fc=i 

The distributions of and are determined by the process {Yk}- We take for {Yfc, fc > 
1} a Markov chain on the two states {F, A} with initial distribution v — {vf,va) and 
transition matrix 

1 — a a 
b 1 - b ' 

where we assume < a, b < 1. The distribution of if„ is then well known, and is called a 
Markov binomial distribution (MBD) (see, e.g., [51 111)). 

Clearly the stationary distribution (np, 7Ta) of the Markov chain {Y^, k > 1} is given by 

b a 
a + o a + o 

It is useful to consider the excentricities e-p and £a of an initial distribution v given by 

v T 

E T :=e T (v) = l , for Tg{F,A}. 

We can then write for k > 1 

Pv(n =r) = 7 r r (l-e r7 fe - 1 ), 
where 7 = 1 — a — 6 is the smallest eigenvalue of P (see also [5J for the computations) . 

4. Probability generating functions of K n and 2f£ 

We compute in this section the probability generating functions of K n and K„- These 
are useful when we consider the convergence of the random variables S n (t) and S^(t) as 
n goes to infinity. 

Given n > 1, let /„ be the probability mass function of K n , i.e., 

f n (j)=F v (K n =j). 

In particular f n (J) = if j < or j > n. Straightforward computations as in [16] or [5] 
yield that 

f n+2 (j + 1) = (1 - 6)/„+i(j + 1) + (1 - a)/„+i(j) - (1 - a - 6)/„(i) 
with initial conditions 

/i(o) = i>a, /i(i) = kp; 

/ 2 (0) = i/a(1-6), / 2 (l) = ^A& + ^a, /a(2) = i*(l-a). 



(4) 



P(F,F) 
P(A,F) 



P(F,A) 
P(A,A) 



Let G n be the probability generating function of K n , i.e., 

71 

G n (s)=E v [s K "] =J2fnU)s j . 

3=0 

It follows from the above recursion equation for /„ that 

G n+2 (s) = ((1 - a)a + (1 - b))G n+1 (a) - (1 - a - b)sG n (s) 

with initial conditions 

Gi(s) — va + vfs, G 2 (s) = va(1 — b) + {vk_b + vva)s + z/p(l — a)s 2 . 

By solving the difference equation of G n with the initial conditions we obtain the proba- 
bility generating function of K n (see also [IS]). 



(5) 



where 



u A (l - P(s) + b(s - 1)) + is F s(a - pjs) + s(l - a)) 

Gn{s} ~ a( S )-m a[s) 

^a(1 — a(s) + b(s — 1)) + v F s(a — a(s) + s(l — a)) _ x 



j8(*)-a(«) 



(6) 



a(s) = i ((1 - a)s + (1 - 6) + \/((l - a).s - (1 - 6)) 2 + 4a&s) , 
= ~ ((1 - <*)s + (1 - 6) - - «)* - (1 - & )) 2 + ■ 



Next we are going to consider the probability generating function of for r G {F, A}. 
Given n > 1, let f£ be the probability mass function of i.e., 

(7) = P* = j) = P* (K n =j\Y n = r). 

In order to deal with it is simpler to deal with the partial probability mass functions 

fa(j) = P„ (K n = j, Y n = r) = fl(j)V v (Y n = t) , 

since these satisfy the same recursion equation as /„: 

/n +2 C? + 1) = (1 - b)f T n+l (j + 1) + (1 - a)/; +1 (i) -(I- a- b)fl(j). 

Only the initial conditions are different: 

A F (0) = 0, A F (l) = ^ F ; 

/ 2 F (0) = 0, fl(l)=v K b, / 2 F (2) = ^ F (l- a ); 

and 

/ a (0) = ^a, A A (1) = 0; 

A A (o)-// A (i-&), #(i) = i*a, A A (2) = o. 
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Then using the above recursion equation of with these initial conditions, the probability 
generating function G T n of can be obtained in a similar way as for G„ (see also [T6]). 

p y (y n = F) tt f (] 



^( s ) = E/n F 0V = E 

(8) *=° 



1 - 

1-EF7"- 1 ) E^ S ' 



j=0 



^ vfP{s) + ^f(1 — a)s 



(a(s)-/3(s)V F (l-£ F7 ™- 1 ) 



TT sa(s) T 



and 



/„ A (J) 



( / 3( s )-a( S )) 7 r F (l- £F7 ™- 1 ) SP[S) ' 



(9) i=o i=o " v " ; j=o 

_ ^a(1 - 6) + ;/ F as - ^Ajg(£) , ^„_i ^a(1 - b) + ^ F as - ^a(s) „, , T 
" (a(s) - /3( S ))vr A (l - eat"" 1 ) " W " «(s)K(l - e Al "-i) PW 



5. Towards continuous time 

To get closer to the PDE model in Section [U we have to fix t = nAt > and then let the 
time step At tend to 0. Hence n goes to infinity. We consider the rates of changes A and 
[i from Section [21 Since the probability that a particle changes its state is proportional to 
the length of the time step At (if At is small) , we should put 



(10) 



, . Xt , at 
XAt = — , b = ixAt = — 



n n 

in the transition matrix P in (|4]). Under this assumption, we will show in this section 
that the random variables S n (t) and S£ (t) defined in Section [3] converge in distribution 
to some random variables S(t) and S T (t) respectively. To achieve this, we first consider 
the characteristic function ip t ,n of S n (t), i.e., 



<Pt,n(u) = E, y 



JuS n (t) 



= E U 



where G n is the generating function of K n given in (JSJ). It is well known that (cf. [2]) 



E^E4(X 1+ ,A^] 



fc=0 



IXi+vAtf] =o(At)=o(-), 
1 J n 



where the last equalities hold since E„ IXf] — o(At) in Equation Q and t is always 
assumed to be fixed. We then obtain by © that 

„,2 



E, 



(11) 

defining 



Ju(X 1 +vAt) 



= 1 + ivE v [X x + vAt] - — E„ [(Xi + vAtf 



n 



tutiv — Du) , 1 s 

— -+° -) 

n n 



z A, 
-+o(- 

n n 



tu(iv — Du). 
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Substituting (fl"0|) and ([TT|) into Equation (j6|) and letting n go to infinity, we obtain 
/(l - + £) + (!-£)+ V((l - + f ) - (1 - £)) + + f ) , 1 A 



+ o - 



exp 
exp 



(i [* - (A + A*)t + ^/(* + t(/i - A)) 2 + 4A/it 2 j ) 

( - 2 [ Dy2 -ivu + X+ Ll- \/{Du 2 -ivu + A - /i) 2 + 4A/jJ) . 



Here we chose the complex square root of (Du 2 — ivu + A — /i) 2 + 4Xu with positive real 
part. 

Similarly, the corresponding limit for /3(E„ ^e m ( Xl +' uAt )] )™ is obtained by replacing the 
minus in front of the square root of the last equality by a plus. It seems convenient to 
introduce the following two notations: 



(12) 



9a = A (u) :=\J (Du 2 -ivu + X- /x) 2 + 4A/i - (Du 2 - ivu), 
9 F = 9 F (u) :=\/(Du 2 - ivu + X- a) 2 + 4A/i + (Du 2 - ivu). 

Then the limits for a(E u [c lu ( Xl+vAt ^] )™ and p(E v [e iu ( Xl+vAt *>] )™ can be rewritten as: 



lim a(E„ 

n— >oo 

lim /3(E„ 



Ju(X 1 +vAt) 



)™=exp((0 A -A-/i)i/2), 
)"=cxp(-(0 F + A + M )</2). 



So we obtain by substituting (fTU|) and (|TT|) into Equation ([3]) that the limit of the char- 
acteristic functions ift >n of S n (t) is a function </? t given by 



(13) Vt (u) 



= (e A -A- M )t/2 ^a6*f + 6> A + A + fj, (eF+x+fl)t/2 v a Oa + v f 9 f -X- h 
9a + f 9 a + 9 f 



It is easy to see that ip t is continuous at u = 0. This implies that there exists a random 
variable, which we call S(t), such that as n — > oo 

S n (t) — > S(t) in distribution. 

Next we are going to consider the convergence of the random variable Sf t (t) as n goes to 
infinity. In a similar way as for S n (t) we consider the characteristic function ipj n of S„(t), 
i.e., 



CLlEy 



jitifXi+oAt) 



where is the probability generating function of K% given in Equation flS)). Substituting 
(|10p and ([TT]) into Equation © and letting n go to infinity, we obtain that the limit of 
the characteristic function cpf n of Sf t (t) is a function ipf given by 

("Ml cn F(,A = r (9 A -A-M)t/2 (^A - X - u) + 2li ( g F+x+fl)t/2 V F (9 F + A + u) - 2u 

1 i4j ^ W (0 A + F Wl-e F .4) +e (0 a + f W1- £f .A)' 
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where A — A(t) — exp (— (A + fi)t). Here we point out that the stationary distribution 
(7tf,7ta) and the excentricities £f,£a do not depend on the time step At, since by fj 10|) 

b fi a A A + /J A + /x 

ttf = — — = — — , 7T A = — — = — — , e F = 1 Up, e A = 1 — v A . 

a + b A + fi a + b a + fi fi A 

Again there exists a random variable, which we call S F (t), such that as n — > oo 

Sl(t) — > S F {t) in distribution. 

Similarly, substituting (|TT))) and (jTTJ) into Equation © and letting n go to infinity, one can 
show that there exists a random variable S A (t) such that S A (t) — > S A (t) in distribution 
as n — » oo, where S A (t) has characteristic function: 

(15) u> A (u) - e (QA-A-M)t/2 ^a(^f - A - /x) + 2A (0>+A+ft)f/2 ^a(^a + A + m ) - 2A 
1 j ^ tlj (0A + e F )irA(l-e A A) (e A + e F )irA(l-e A A)- 

6. Modeling the kinetics with a continuous time Markov chain 

In our model we used a simple discrete time set up. This will be useful in Section [9l 
but it is worthwhile to compare our results with a model that involves a continuous time 
Markov chain. Let Y(t),t > denote the state of the particle at time t. Recall from 
Section that A and fi are the rates of changes from 'free' to 'adsorbed' and 'adsorbed' 
to 'free' respectively. Hence it is natural to model the kinetics by a two-state continuous 
time Markov chain {Y(t),t > 0} with initial distribution P„ (F(0) = r) = v Ti t e {F, A} 
and generator matrix 



Q = 



Q(F,F) Q(F,A) \ = / -A A 
Q(A,F) Q(A, A) ) \ fi -ii 



The solute can only move when it is free, and in this case we model the displacement due 
to dispersion and advection as a Brownian motion with drift v. 

A trick to deal with continuous time Markov chains is uniformization. This idea gives us 
an alternative way to model the S(t) and S T (t),T £ {F,A} obtained in Section [5j Let 
A > max (A, fi) be the rate of the uniformization. It follows that (see e.g. [13], page 402) 
the continuous time Markov chain {Y(t), t > 0} can be viewed as a discrete time Markov 
chain {Zk,k > 0} over the same state space {F,A} and the same initial distribution 
P u (Zq = t) = v T , t £ {F, A}, but with the transition matrix 



Pa = 




P(F,A) 
P(A,A) 




Let N(t) be the number of the state transitions up to time t, which is a Poisson process 
with rate A. Let ifjv(t) be the occupation time of the chain {Z^} in state F up to time 
N(t), which is a Markov binomial distributed random variable, when conditioned on N(t). 
Since the solute can only move when it is free and the displacement in the free state is due 
to dispersion and advection, we model Afc, the displacement during the fcth free interval, 



10 



as a Brownian motion with drift v stopped at time T which is exponentially A distributed. 
So we put 

X k = Af(vT, 2DT) with T = Exp(A). 

Then we can write T-L A {t), the position of the particle at time t with respect to the 
uniformization at rate A, as: 

fc=l 



Similarly, for r G {F, A} we can define H A (t) merely by changing -Kjv(t) to the conditional 
random variable KJfU) = ^iv(t) I {%N{t) — r }i i- e -i %a(*) denotes the position of the 
particle at time t conditioned on being in state r at time N(t). 

Letting A go to infinity, one can show by using the characteristic functions of Ti. A (t) and 
T-L A (t) as m Section [5] that 



H A (t) — > S(t), n T A {t) — > S T (t) in distribution, 



where S(t), S T (t) are the same random variables as in Section [SJ 

It is even more natural to look at the continuous time Markov chain {Y(i), t > 0} directly. 
Let U (t) be the occupation time of the chain {Y(t)} in state F up to time t, and let fu(t) be 
its probability density function. We model the displacement of the solute in the free phase 
as a Brownian motion with drift v. Then the position H (t) of the particle at time t can 
be written as a normal distribution with mean vU(t) and variance 2DU(t). Conditional 
on U (t) it follows from Equation (5) of [12] and Equation ((T3|) that for <f> > 



■iH(t) 



OO pOO 



o Jo 



^—(Du —ivu)x —<fit 



fu(t) (x) dx dt 



A + (j) + Du 2 — ivu —A 



4> + A + + v\(Du 2 — ivu) 
((f> + \ + Du 2 — ivu)((f) + /i) — \/i J v 



a iuS(t) 



-^dt. 



Since E^ [e™ 5 '- 4 '] is a continuous function of i by Equation (jT5)l . it follows from Lerch's 
theorem (cf. [14], page 24) that E^[e raff W] = E„[e m5(t) ] for all t > 0. Hence H(t) and 
S'(t) have the same distribution. 

Similarly, for t e {F, A} let H T (t) = H(t) \ {Y(t) = r} be the conditional random variable 
denoting the position of the particle at time t conditioned on being in state r at time t. 
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From the proof of Theorem 1 in [3] one obtains that for <j> > 
P„ (Y(t) = F)E U 



~ e iuH F (t)~ 


pOO 

e^*di = / E v 




Jo 



e y> l {Y (t)=F} 



A + <f> + Du 2 — ivu 



fJL + ■ 



A + Du 2 — ivu){<j) + fi) — A/i 



P„ (Y(t) = F) E„ 



iS F (t) 



-^dt, 



where the last equality follows using Equation (|14p and since the limiting probability of 
a particle being in state r at time t is given by 

(16) P„ {Y{t) =t)= lim P„ (F„ = r) = lim vr T (l - £ T 7 n_1 ) = tt t (1 - e T A(t)) 

n— >oo n— ^oo 

with A(t) — exp (—(A + n)t). Similarly one can also show that for <fr > 



p„ (y(t) = A) e„ 



-^dt 



A + ^a(0 + -Dit 2 — tra) 
+ A + Du 2 — ivu){4> + A*) — A/i 



p„ (y(t) = A) e„ 



juS A {t) 



-^dt. 



Again, using Lerch's theorem, it follows that H T (t) and S T (t) have the same distribution. 
Therefore our discrete time model converges in distribution to the same random variables 
as obtained by the natural continuous time Markov chain. 



7. Densities and partial differential equations 

We will show in this section that for instantaneous injection of the solute, i.e., with initial 
distribution v — (1,0), the partial probability density functions fg(t,x) and fg(t,x) of 
S F (t) and S A (t) do satisfy the partial differential equations in (jTJ) . 

Let fs{t,x) and fg(t,x) denote respectively the probability density functions of S(t) and 
S T (t) for r 6 {F, A}. Recall from (|16l) that the probability of a particle being in state r 
at time t is given by P^ (Y(t) = r) = 7iy(l — e T A(t)). We define the partial probability 
density functions of S T (t) as 

(17) f T s {t,x) = P„ (Y(t) = r) r s {t,x) = tt t (1 - e r ^(t))/^(i,x). 

Obviously f s (t,x) = f*(t,x) + f£(t,x). 

Lemma 7.1. Lei 0& = &a(u), Of = #f(m) 6e defined as in (IT2|) . TTien 

lim |u 2 (6>a(w) - (A - //)) I = 2Xfi/D, lim |6> F (w)/w 2 | = 2D. 
Proof. It is straightforward to check these formulas. □ 
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Lemma 7.2. The probability density function fg(t,-) of S F (t) can be written 

f s {t,x) = ±- J e-™rf(u)du, 
where tpf is the characteristic function of S F (t) given in (|14p . 



Proof. We only need to show that <p F is integrable. Obviously <f F is a continuous function. 
So it suffices to show that Ji u i >m \<ft du < oo for some M > 0. From Lemma ITTTl and 
(|14[) it follows that for all |m| large 



(0 A + ^FhF(l-e F ^) 



(e F +A+ At )t/2 M^f + A + /x) - 2/x 
(^A + ^F)7r F (l-£F^) 



<^ + C 2 e-^ 2 / 2 , 



□ 



where C%, C 2 are constants independent of u. This finishes the proof of the lemma. 
Surprisingly, Lemma T7.2I does not hold for S A (t), but we still have the following. 
Lemma 7.3. The distribution [i A of the random variable S A (t) can be written as 

[i A = kS + (1 - n)fi A 

where n — v A e~ ^ / {-k A (l — e A A)) and jl A is the distribution of a continuous random 
variable having probability density function 

1 



m,x) 



2tt(1 - k) 

with (p A the characteristic function of S A (t) defined in U5)) . 

Proof. It follows from Lemma \7. II and Equation (fT5|) that for all |u| large 



<pH u ) - K 



< 



a (e A -A-^)t/2 ^a(6'f-A-^) + 2A 



v A e 



-fit 



(18) 



< 



(0 A + F )7r A (l-£ A .A) n A (l-e A A) 
-(g F +A+M)t/2 ^a(Aa + A + ^)-2A 
{0 A + e F )ir A (l-e A A) 
je A -x-u.)t/2 2A-t/ A (^A + A + /z) 



v A e 



- fJL I 



7T A (1 - e A A) ^ 



+ 6» F )7r A (l-eA^) 



,(0 A -A+/i)t/2 



+ C 2 c- Dtu2 / 2 



where C\ , C 2 are constants independent of u. This implies that the integrand in the lemma 
is integrable. 

Without loss of generality we may suppose v A > 0. Using (fl8|) we obtain that asT->oo 



2^ f T l P J t( u ) du - K = J (ft ( U )~ K ) du 



0. 
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This implies that the point is an atom of /l«a, since (cf. 2], page 306) 



Ma({0}) 



lim — 

T^oo 2T 



<Pt (u)Au = K. 



Moreover, is the unique atom of //a since (cf. [2J, page 306) 

5>a({?})) 2 = lim ^ / T |^H| 2 d M = ( MA ({0})) 2 , 



where the sum is taken over the set of points of positive /xa measure, and the second 
equality can be seen by using (fT8]l and the fact that <fit~(u) is uniformly bounded. This 
establishes the lemma. □ 

It follows from Lemma 17.31 that S A (t) is a continuous random variable if and only if 
v = (1, 0), i.e, for instantaneous injection of the solute. It is interesting that in this case 
we have the following. 

Theorem 7.1. The partial probability density functions fg of S T (t) for t € {F, A} satisfy 
the partial differential equations {1} : 

dfs&x) , 9fj(t,x) _ D d 2 fl{t,x) JfE&x) 



dt 



dt 
dt 



dx 2 dx 
-fifi(t,x)+\fUt,x) 



for t > 0, with initial and boundary conditions 
fi(0,x):=S(x), /#(0,x):=0; 

lim f s (t,x)= lim d -M^- 



for t > 0, t G {F, A}. 



Proof. The initial conditions imply v — (1,0). It follows from Lemma T7. 21 17.31 and Equa- 
tion {32} that 



e- mx ifl{u)du for re{F,A}, 



where 
(19) 



<p T t {u) = 7T r (l - e T A{t))ip T t {u), 



with ipj the characteristic functions of S T (t) given in (fl4|) and (|15|) respectively. It is easy 
to see that fg and fg satisfy the initial and boundary conditions. 
Using Lemma 1 7. II it is not hard to check that the four functions in u 



de lux (p^{u) 




de- iux <p£(u) 




de tux (p^(u) 




d 2 e- iux $(u) 


dt 




dt 


7 


dx 




dx 2 



are all bounded by a function of the form C\ju 2 + C2C 2 Du for \u\ large, where C\, Ci are 
constants independent of u. Thus we can exchange the integral and differential operators 
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in the partial differential equations (cf. [6], page 417). Hence we only need to show that 



;A 



-Z?u t (u) + ivutpl (u) — \ip t (u) + [i(p t (u) 



^M = _^(„)+A^(«). 
It follows from (JHJ), (JTSJ), O and ^ A = that 

— Du 2 ipf(u) + ivuipj(u) — Xipf (u) + fiipf{u) 

_ 0A-&F _ (g A -A-g)t/2 gA ~ A + ^ | p -(9 F +A+M)t/2 g F + A-/i \ 



2A^ 



A + #F 

_ fl A - A - H c (0 A - X -n)t/2 °A - A + /i _ 6> F + A + M c -(fl F + A +M )t/2 6> F + A - /X 

2 6>a + 6>f 2 1 6>a + 6» f 

St ' 

where the second equality holds since (9a — A + /i)(# F + A — fj,) = 4A/x. 
One finishes the proof of the theorem by checking 

d JfM = -^ { u) + X^{u). 

□ 

We would like to point out that Lindstrom and Narasimhan .9 gave an analytical solution 
of the partial differential equations with different initial and boundary conditions by using 
Laplace and inverse Laplace transforms. Their method can also be used with our initial 
and boundary conditions to give the same solutions as we have obtained via our stochastic 
model as Theorem 1 7. II 



8. Moments of S(t) and S T (t) 

The mean and variance of S(t) can be obtained by differentiating its characteristic function 
ipt given in (|13p . but a more leisurely way is to take the limits olE v [S n (t)] and Vax v (S n (t)) 
respectively. 

Lemma 8.1. The first and second moments of S(t) can be obtained by taking the limits 
of the corresponding moments of S n (t) respectively, i.e., 

E„[S»(t)] -> E„[S(t)] , Var„(S n (t)) -► Var„(S(t)) 

Proof. Recall that the mean of K n is given in [5] . It is not difficult to check that the first 
moment of K n is uniformly bounded, i.e., there exists M > 0, such that |E v [if n ] | < M. 
Since the Xk's are independent random variables also independent of K n , using At — t/n 
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and ([3]) we obtain that 

71 

E„ [S 3 n ] = E„ [E„ [S 3 | iC„]] = /"0') E " 

j=o 



(^(X fc +«At))' 



= E„ [(Xi + ^At) 3 ] £ /„(j)j + 3E„ [(Xi + vAt) 2 } Eu[X! + vAt] ^ f n (j)j(j - 1) 

2=0 3=1 
n 

+ (E„[JTi + vAt] ) 3 g - 1)0" - 2) 

< (t + 6D^ 2 + « 3 t 3 )E„[7f n ] + 3(2L>t + w 2 t 2 )wtE„[.K„] + v 3 t 3 , 

which implies that S n (t) and S^(t) are uniformly integrable. This together with the fact 
that S n (t) converges to S(t) in distribution (shown in Section [5} imply that E„ [£„(*)] — > 
E v [S(t)],Var v (S n (t)) ->■ Vaiv(S(t)) (see e.g. 0, Theorem 25.12). □ 

Since Xk is independent of K n , from Equation (4) in [3] and Proposition 2.1 in [5] together 
with Equation (J3]) we can determine the first and second moments of S n (t): 
E v [S n (t)\ = E v [K n ] E v [Xi + vAt] = E v [K n ] vAt 

( e F (l- 7 "A £F 7^(1-7") 

= TTp [ 71 VAt = 7TF v t vAt, 

V 1 -7 / 1-7 

and 

Var v (S n (t)) = E„[i<r n ] Var^ + vAt) + Vax u {K n ) (E U [X 1 + vAt}) 2 



7T F 



( n _ EEkL 2 DAi + Var„(i<f n ) (wAi) 2 



= 2D** i - 2^(1- 7") Af + -A(l + 7) + 2 £F(7 r A -, F ) 7 ^ F ^ 
1 — 7 1 — 7 

7 (eF(7rF-7rA)-27r A ) -e F (7TA-^F) , A ,, 2 

+ ^—^y 2 7T F («A<) 

„/£ F (7r F -tta) 7tta + e F (7r A - J^f) n 7r F e 2 \ 2 

Substituting At — t/n and (|10[) in the mean and variance of S n (t) and letting n — > oo, by 
Lemma 18. II we obtain the mean and variance of S(t). 

Proposition 8.1. The mean and variance of S(t) are given by 



and 



E v [S(t)]=w F vt-^-v(l-A), 
A + fi 



,n,.^ „^ , 2De F 7r F , 2(7r A + e F (7r A - n F )A) 

Vax v (S(t)) = 2Dw F t (1 -A) + — r- -^-ir F v 2 t 

A + /i A + fi 



e F (7r F - tt a ) - 2tt a - e F (7r A - v v ) 



-n-pv 



2 



(A + M) 2 

V (A + a*) 2 (A + [i) 2 J 



2 
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Now we are going to consider the means and variances of S T (t),T € {F, A}. Again, one 
could obtain them from their characteristic functions, but we will use the following lemma, 
which can be proved in a similar way as Lemma 18.11 

Lemma 8.2. The first and second moments of S T (t) , r £ {F,A} can be obtained by taking 
the limits of the corresponding moments of S^(t), i.e., 

E v [S£{t)] -> K[S T (t)} , Vfa v ffi(t)) -> V&T v (S T (t)) . 

Because of independence, using Equation (5) in [5] and Proposition 3.1 in [5] together with 
Equation ([3]), we obtain that 

P rcPml v Tr-fia, TTp-gFTTAT"" 1 , , (7r A -£F7r F )(l-7") A , 

and 

VeXv(Sl(t)) = E„[<] Var„(X fc + vAt) + Var„«) (E„[X fe + vAt]f 
1-£f7™~ 1 (1 - 7)(1 - £f7" _1 ) 1-EF7™- 1 

7T F -£ F 7TA7"~ 1 £ , (7TA -£ F 7T F )(1 ~7") A A ^2 

l-£ F 7 n_1 (1 - 7 )(1 - £ F 7 rl_1 ) 

7r A 7r F (l + 3£F7 n ~ 1 ) , 2 £ F ^+7rl7"-2^ A ^ F (l + £ F7 "- 1 ) \ 



1-EP7™" 1 (1 - 7 )(1 - EP7"- 1 ) 

, ^ 7r A 7r F (4 + £ F ) - (tt a + £ F 7rg) e F tt| + tt a - 2tt a tt f (1 + eg) \, ,.x 2 
+ (1 - 7 } { (1 - 7)(1 - £F7- 1 ) + 2 (1 - 7) 2 (1 - £F7- 1 ) ) (vAt) ■ 
Substituting At = t/n and (|10p in the mean and variance of S F (t) and letting n — ► oo, by 
Lemma 18.21 we obtain the mean and variance of S F (t) . 

Proposition 8.2. The mean and variance of S F (t) are given by 



_ 7Tf — £ F 7TX«A_ iit , (TTA - £ F 7Tf)(1 - -4) 

(A + m)(1-e f ^) 



E„ \S* (t)\ = — Vt + ^ 7-7^ rr+V, 

1 — E F J\ 



Y^(S F (t)) = 7TF - £F7TA A A 2Dt+ - *A- e Ftt* ^(i, ,4) + *f " 

v y l-e F -4 (A + ^)(l-£ F -4.) 1-ef-4 



7T F -£ F 7TA^l, , 7T A - £ F 7T F \ 2 2 

X + u) 1 - s F A) y ' 



l-e ¥ A (X + fi)(l-e F A) 
n e F nj + n%A - 2tt a tt f (1 + e F A) 2 . _ £ F tt f + tt a - 27r A 7r F (l + e F ) 2 
(A + /x)(l-£ F .4) + ^ (A + M ) 2 (1-£f-4) 

In a quite similar way we obtain the following result. 

Proposition 8.3. The mean and variance of S A (t) are given by 

[sA(] = n F -e A n A A v (e A n A -n F )(l-A) 
1 WJ l-s A A (A + /i)(l - e A A) 
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and 

Var,(5 A (t)) = **- e ^ A 2Dt + ^a - tt f _ tt| - ea^M 2f2 

v y 1-eaA (A + /i)(l -ea-4.) I-ea-4 

7T F - EA^aA £ A 7T A - 7TF ,.. _ ,,\ 2 2 

(A + M )(1-£ A V " 
^ Trg + £ A 7T A .4 - 7r A 7r F (l + £ A )(1 + A) 2. oh _ ,, 7rg +£ A 7T A - 27r A 7T F (l +£ A ) 2 

(A + A*)(l- eA ^) [ ] (\ + n) 2 (l-e A A) V - 

Now, we will use our model to illustrate some mistake made by Michalak and Kitanidis 
|10j . Note that the moments of S F (t) and S A (t) calculated in Proposition 18.21 and 18.3 
are exactly the same as that calculated by Michalak and Kitanidis (see Section [2]) if the 
initial conditions are specified as v = (f,0), since for m > 1 

E ( i,c»[(^(i)n = / x m m,x)dx = [ X ™^A& X = M<rHt\ T & {f, a} 



where the second equality holds since by Theorem 1 7 . 1 1 both the partial probability density 
functions fg = 7r T (l — e T A)fg and the concentration functions C T satisfy the partial 
differential equations (TT|) with the same initial and boundary conditions. 
Recall from Section [2] that we translate the parameters into our paper as follows: 

jji = k, A = /3k. 

If we let the solute be 'free ' at time and t, i.e., the initial distribution v — (1, 0), then 

A Ml A (3 

£f = = -p, t^f = t— — = a ■ - , vr A = 



H X + fi + r A + /z (3 + 1 

Substituting these parameters into Proposition 18.21 yields 



1 + 0* A 2,2 / 1+^4 2/3(1 - ^l) ^ 

V t - [ — —- — — t + — — —TT— — - V 



G8 + l) 2 (l + M) \{/3 + l){l + f3A) k(p + l) 2 (l + /3A) 

6(3((3A~l) „ 6008-1) 2 , 



.2 



k(J3 + 1) 3 (1 + (3A) V t+ k 2 (/3 + 1) 4 (1 + PA) V (1 ^' 

where .4 = cxp (—(A + fj,)t) = exp (— (/3 + l)fct). This gives indeed Equation © which is 
taken from |10) . 

However, Michalak and Kitanidis state in their paper that Va.r„(S T (t)) can be obtained 
by a linear combination of Var p (S' T (t)) and Var A (S ,r (t)) (i.e., Vax v (S T (t)) with initial 
distributions f = (1,0) and v — (0,1)). This is not true, and we provide the correct 
formulas in Proposition 18.21 and 18.31 We have also given the formula for the total solute 
in Proposition 18. II 
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9. Double-peak behavior in reactive transport models 

Double peaks in the 'free' concentration distribution Cf are discussed by Michalak and 
Kitanidis [TU] using simulations. Theorem 17.11 tells us that Cf(£, •) can be seen as the 
partial probability density function fg(t, •) of S F (t) if the initial distribution is v — (1,0). 
We will show in this section how double peaks can also be explained by means of our 
stochastic reactive transport model. Let f$ (t, •) be the probability density function of 
S F (t) defined in Section[3] We are going to approximate fg (t, ■) by fg (t, ■), since S F (t) 
converges to S F (t) in distribution. 

Michalak and Kitanidis consider Gaussian diffusion, i.e., the Xk's are normally distributed 
random variables with mean and variance 2D At, which satisfy Equation ([3]). So the 
characteristic function of S F (t) can be written as 

n 

) = fnU) exp (ivAtju - DAtju 2 ) , 

3=0 

where f F is the probability mass function of K F defined in Equation ([7]). Obviously 
I^°co \*^t n( u )\du < oo. Thus by the inverse Fourier transformation, using that f F (0) = 0, 
we obtain 

1 r°° 

/ e- ia *<pl n (u)du 

n ^ „oo 

(20) = Z2 /nW^ J ex P [HjvAt -x)- u 2 jDAtjdu 

V^,f,-n 1 ( (x~jvAt) 2 \ 

= L fn ^^JBAl eXP ( AjDAt ) ■ 

So Sl(t) is a mixture of Gaussian distributions with mean jvAt and variance 2jDAt. 
Recall from [5] that the probability mass function f F of can be unimodal or bimodal. 
This property of K F gives rise to the same phenomenon for Sl(t), i.e., one peak or two 
peaks appear in the probability density function fg n (t, x) of S F = S F (t) for large n. 
Michalak and Kitanidis focus on the case that the solute starts in the free phase and the 
length of the initial solute is L, i.e., the initial distributions of the PDE's (|T|) are given by 

C F (0,x) = ^l [0tL] (x), C A (0 1 x) = 0. 

So to make the comparison, we look at the probability density function /| (t, x) of 

sl(t) = sl(t) + u L , 

where Ul is a uniformly distributed random variable over the interval [0, L] (independent 
of S F (t)). Michalak and Kitanidis point out that the double peaking behavior of the free 
concentration distribution is a function of the so called Damkohler number of the first kind 
Dai — fxLR/v, where R is the dimensionless retardation coefficient. They state that the 
timing of its appearance is controlled by the mass transfer rate and the retardation factor, 



<Pt,n( u ) = E » 



JuS*(t) 



= G 
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FIGURE 3. The three graphs in the left column are the normalized con- 
centration functions Cp(t, -)/max x C?F(t,x) copied from Michalak and 
Kitanidis [10 . The three graphs in the right column are the normal- 
ized probability density functions /| (t, -)/max x f^ (t,x) given by the 
Fourier transformation in our paper. All graphs have Pe = 100, v = 
L = 1, R = 2. In the first row Daj = 0.1, t* = 3.6, in the second row 
Da T = 0.33, t* = 3.2, and in the last row Da/ = 1.0, t* = 3.0. 



i.e., the dimensionless time t* = [i(R — l)t. The so called Peclet number Pe = vL/D is 
kept constant at a value of 100. Recalling from Section [2] that A = fin = (R — l)/x and 
a = XAt, b = fj,At, we translate these parameters into our paper as follows: 

k(R-l)t t* , kt t* vL . t t*LR 

a=— '- = —, b= — = — -, D = —, At=- = — — — — . 

n n n n{R—l) Pe n nv(R—l)Dai 

The graphs in the left column of Figure [3] are a copy of the graphs of the normalized 
aqueous concentration functions Cp(t, -)/max x C-p(t, x) (consisting of the free particles) 
in Michalak and Kitanidis [10; using simulations corresponding to different choices of the 
Damkohler number Dai and dimensionless time t* . The three graphs in the right column 
of Figure [3] are the normalized density functions f~ (£,•)/ max x /| (t, x) calculated 

t>400 i>400 

using Equation ([20)) corresponding to the same choice of Dai and t*. The number n is 
chosen large enough such that max (a, d) = max (A At, /x At) < 0.01. From Figure [3] it is 
obvious that our model gives a much better view at the double peaking phenomenon. 
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/ \ Da, = 0.46 




/ / V \ Da, =0.3 






Da, =0.1 

^ ~~ - ^ / \ 



20 40 60 80 

FIGURE 4. The three graphs are the probability density functions 
/f (V) of Sj 00 (t). All graphs have Pe = 100, v = L = 1,R = 2,t* = 
4.0, and different Damkohler numbers. 



Moreover, for each t* , by a numerical calculation we can obtain upper bounds for Daj 
such that double peaks appear. For example, Figure HI gives an intuition on how double 
peaks behave when Dai increases. We numerically calculated the upper bounds for Dai 
in Table 1 corresponding to different dimensionless times t* with R — 2. For example, 
when t* = 2.0 two peaks occur for all Dai > until Dai = Daf ax = 0.43. Table 1 
suggests that double peaking is pronounced for 2 < t* < 5, and almost dies out when 
t* < 1 .5 or t* > 10. 

Table 1. R = 2,Pe = 100, v = 1, L = 1, n = 400. 



t* 


1.5 


2.0 


2.5 


3.0 


3.5 


4.0 


4.5 


5.0 


6.0 


7.0 


8.0 


9.0 


10.0 




0.12 


0.43 


1.45 


1.42 


0.73 


0.45 


0.30 


0.21 


0.11 


0.07 


0.04 


0.02 


0.02 



10. Final remarks 

We emphasize that the so called 'random walk method' or 'particle tracking method' first 
proposed by Kinzelbach [8] has a relation to our model, but has always been used as a 
simulation tool, to perform numerical experiments (for a recent example see pQ). In fact 
it is shown in |15j for the first time that if one takes an appropriate limit (in a similar way 
as in [1]), then the Fokker-Planck equations of an extended version of our simple model 
to a Markov chain which also involves discrete steps in space, yield the partial differential 
equations ([I} in Section [5J 

Finally we mention that our computations yield the following. If one starts in the sta- 
tionary distribution, i.e., v = (7tf,7ta), then ef = £a — 0. Substituting 



e F = 0, 



7T F 



A + ^i' 



7TA 



in Proposition 18. 11 we obtain 
Var„(S(t)) =^^t + 



2fi\ 



X + H (A + fi) 



rvH- 



2fiX 



A + ^i 



(A + /i) 4 



! (l-exp(-(A + A1 )t)). 
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We then recuperate a (more general and more detailed) version of the main result of Gut 
and Ahlberg ([7 , p.251). 
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